!dk wates
      subroutine watec(xp, yp, zp, w) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      real(double) , intent(in) :: xp 
      real(double) , intent(in) :: yp 
      real(double) , intent(in) :: zp 
      real(double) , intent(inout) :: w(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: i, j, k, l 
      real(double) :: nu, the, zeta, wi, wim, wip, wj, wjm, wjp, wk, wkm, wkp, &
         wtotl 
!-----------------------------------------------
!
      i = int(xp) 
      j = int(yp) 
      k = int(zp) 
!
      the = xp - i - 0.5 
      zeta = yp - j - 0.5 
      nu = zp - k - 0.5 
!
      wi = 0.75 - the**2 
      wim = 0.5*(0.5 - the)**2 
      wip = 0.5*(0.5 + the)**2 
!
      wj = 0.75 - zeta**2 
      wjm = 0.5*(0.5 - zeta)**2 
      wjp = 0.5*(0.5 + zeta)**2 
!
      wk = 0.75 - nu**2 
      wkm = 0.5*(0.5 - nu)**2 
      wkp = 0.5*(0.5 + nu)**2 
!
!     k-plane
!
      w(1) = wi*wj*wk 
      w(2) = wip*wj*wk 
      w(3) = wip*wjp*wk 
      w(4) = wi*wjp*wk 
      w(5) = wim*wjp*wk 
      w(6) = wim*wj*wk 
      w(7) = wim*wjm*wk 
      w(8) = wi*wjm*wk 
      w(9) = wip*wjm*wk 
!
!     k-1 - plane
!
      w(10) = wi*wj*wkm 
      w(11) = wip*wj*wkm 
      w(12) = wip*wjp*wkm 
      w(13) = wi*wjp*wkm 
      w(14) = wim*wjp*wkm 
      w(15) = wim*wj*wkm 
      w(16) = wim*wjm*wkm 
      w(17) = wi*wjm*wkm 
      w(18) = wip*wjm*wkm 
!
!     k+1 - plane
!
      w(19) = wi*wj*wkp 
      w(20) = wip*wj*wkp 
      w(21) = wip*wjp*wkp 
      w(22) = wi*wjp*wkp 
      w(23) = wim*wjp*wkp 
      w(24) = wim*wj*wkp 
      w(25) = wim*wjm*wkp 
      w(26) = wi*wjm*wkp 
      w(27) = wip*wjm*wkp 
!
      wtotl = sum(w(:27)) 
!
      if (abs(wtotl - 1.) > 1.E-10) write (59, 100) wtotl 
  100 format(1p,e12.3) 
      return  
      end subroutine watec 
